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Abstract 

A general quantum many-body theory in configuration space is developed 
by extending the traditional coupled-cluster method (CCM) to a variational 
formalism. Two independent sets (destruction and creation sets) of distribu- 
tion functions are introduced to evaluate the Hamiltonian expectation. An 
algebraic technique for calculating these distribution functions via two self- 
consistent sets of equations is given. By comparing with the traditional CCM 
and with Arponen's extension, it is shown that the former is equivalent to a 
linear approximation to one set of distribution functions and the later is equiv- 
alent to a (generalized) random-phase approximation to it. In additional to 
these two approximations, other higher-order approximation schemes within 
the new formalism are also discussed. As a demonstration, we apply this 

technique to a quantum antiferromagnetic spin model. 
PACS numbers: 31.15.Dv, 75.10.Jm 
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I. INTRODUCTION 



The main task of a microscopic quantum many-body theory is to study correlations 
between the constituent particles of a quantum system in a systematic way. The treatment 
of these many-body correlations is either in real space or in configuration space. A real-space 
theory usually focuses on the potential part of many-body Hamiltonians; a configuration 
space theory often starts from the kinetic part of Hamiltonians. One of the most successful 
real space quantum many-body theories is the method of correlated basis functions (CBF) |IJ 
in which real-space correlation functions of the ground state are determined variationally. 
Perhaps, the closest counterpart of configuration space theories to the real space CBF is 
the coupled-cluster method (CCM) ||-f§] in which correlation operators are employed to 
construct the ground state. One key feature of the CCM is that the bra and ket states are 
not manifestly hermitian to one another ||. 

In this paper we propose a general variational theory in configuration space by extending 
the traditional CCM to a variational formalism in which ket and bra states are hermitian to 
one another. The difficult task of evaluating the Hamiltonian expectation can be done by 
introducing distribution functions which can then be determined either by a diagrammatic 
technique or by an algebraic one. The diagrammatic approach developed in this context 
is quite similar to that of the CBF. In the algebraic approach, one derives two similar 
sets of self-consistent equations for the distribution functions; these equations can then be 
tackled by various methods, e.g. iterative method. Easy comparison can be made with the 
traditional CCM in this approach. We will mainly discuss the algebraic approach in this 
article; the diagrammatic approach will be discussed elsewhere || . We apply this variational 
method to a well-known spin model as a demonstration. Some of our preliminary results 
has been reported in a conference paper |7j. 
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II. THE REPRESENTATION OF A MANY-BODY WAVE FUNCTION 



Similar to the method of correlated basis functions (CBF), the coupled-cluster method 
(CCM) deals directly with the wave functions of a many-body system. We shall take the 
spin-1/2 antiferromagnetic XXZ model on a bipartite lattice as an example. The model 
Hamiltonian is given by 

H = 2 E H i,i+p = \ E (Asfsf +p + \s+si +p + l -sist + \ , (1) 

i,p i,p 

where A is the anisotropy, the index / runs over all lattice sites, p runs over all nearest- 
neighbor sites, and s ± are the usual spin raising (+) and lowering (— ) operators. The 
Hamiltonian of Eq. (1) at A = 1 corresponds to the isotropic Heisenberg model which has 
been a focus of theoretical study in recent years due to its relevance to high-temperature 
superconductivity. 

In the limit A — > oo, the ground state of Eq. (1) is clearly given by the classical Neel 
state with alternating spin-up and spin-down sublattices. We shall exclusively use index i 
for the spin-up sublattice and the index j for the spin-down sublattice. For a finite value 
of A, such as the isotropic point A = 1, the many-spin correlations in its ground state can 
then be included by considering the excited states with respect to the uncorrelated Neel 
model state. These excited states are constructed by applying the so-called configuration 
creation operators C\ to the Neel model state with the nominal index / labelling these 
operators. In our spin model, the operators C\ are given by any combination of the spin-flip 
operators to the Neel state, namely and s~j and the index / in this case corresponds to 
the collection of the lattice indices (i's and j's). The hermitian conjugate operators of C\ 
are the configuration destruction operator Cj, given by any combination of sf and sj . For 
example, the two-spin flip creation operator is given by Cjj = s^Sj~, and their destruction 
counterpart, CV,- = sfsj. 

The traditional CCM is based on the Hubbard, Hugenholtz and Coester representation 
(HHC) for the ground ket state, where the correlations are parametrized by an exponen- 
tiated operator as, 
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|^)=e s |$), S = J>/Cj. (2) 

For our spin model, |$) is the Neel state and Fj are the correlation coefficients. The 
configuration creation operator Cj in this case is given by a product of any number of pairs 
of the spin-flip operators, 

N/2 - - + + 

ew) = e E ./v...,,.. 1 ■ w 

I n=lti...Ji... v ZA / 

where s is the spin quantum number. Although we are mainly interested in s = 1/2, we 
keep the factor of l/2s for the purpose of comparison with the large- s expansion. Notice 
also that in Eq. (3) the spin-flip operators of the 2-sublattice always pair with that of the 
j-sublattice to ensure the total z-component s z total = 0. For the bra state, however, the CCM 
proposes a different, practical form as [0-0, 

(tf | = ($|S'e- s , (4) 

where S is as given in the ket state and the linear bra state operator 5" is constructed by 
the configuration destruction operators only, namely, 

N/2 + + - 

#=i+£^ci = i+£ e iv " s ;;^" S3n • (5) 

I n=lti...,ii... v b > 

The coefficients {Fj,Fj} = {f^...^..., fi 1 ...j 1 ...} are determined variationally through the 
Hamiltonian expectation (H), noticing the normalization condition (^|^) = 1, 

(V\H\V) = (®\H\<$>), (6) 

where the similarity-transformed H = e~ s He s can be expanded as a series of nested com- 
mutators as 

H = H + ±[H,S\ + ± i [[H t S\,S\ + ..- • (7) 

In most cases, H contains a finite order of destruction operators. The above series then 
terminates at a finite order as S contains only the creation operators. Hence, the Hamilto- 
nian expectation value in the CCM is a finite order polynomial function of the coefficients 



{Fi,Fi}. More specifically, (H) in the CCM is linear in the bra-state coefficients Fi and 
finite-order polynomial in the ket-state coefficients Fi. Thus, calculations in the CCM in 
general are quite straightforward; when an approximation scheme is chosen (i.e., a trun- 
cation scheme with a finite set of {Fj,Fj}), no further approximation is necessary in most 
calculations. However, this CCM parametrization of the ground state is problematic in deal- 
ing with long-range correlations as discussed in the context of our spin model calculations 
in Ref. ||. More discussion of the problems in this traditional CCM will be given later. 

An obvious extension of the CCM is to apply the HHC representation to both the ket 
and bra states. Hence we write 

\^) = e s \$), S = Y,FiCl (8) 
(9\ = (*\eP, S = Y,FiCi, (9) 

where the ket-state and bra-state correlation coefficients Fj and Fj are hermitian and inde- 
pendent to one another. For our spin model, the model state |<&) is the Neel state, and the 
correlation operators J2i FjC\ and J2i FjCj are given as in Eqs. (3) and (5). The coefficients 
{Fj, Fj} are then determined by the usual variational equations as 

SFj ori 

where energy expectation is defined in the usual way as 

= (£|g|g) = (gjefgefjg) 
K ' ($|e s e s |$) ' 

Clearly, the normalization factor (^|^) and Hamiltonian expectation (H) are highly non- 
trivial functions of the coefficients Fj and Fj. Their calculation in the standard variational 
approach is in general difficult, contrast to the CCM where the expectation value of the 
Hamiltonian is a finite-order polynomial of the coefficients, as described earlier. This is 
perhaps the main reason that little progress has been made in this extension of the CCM 
except some general discussion and a few attempts M. 
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Hence, the key to the extension of the CCM to the standard variational method as de- 
scribed in Eqs. (8)-(ll) is to develop a practical and consistent technique to evaluate the 
normalization factor and the Hamiltonian expectation. It is known in statistical mechanics 
and real space quantum many-body theory that these evaluations can be done more effi- 
ciently by employing distribution functions. One then needs to develop a systematic and 
consistent scheme to calculate these distribution functions. We have considered two such 
schemes. One is similar to the traditional technique in statistical mechanics employed by the 
CBF. In this method, one introduces a generating functional whose functional derivatives 
are the distribution functions. A diagrammatic technique has been developed to evaluate 
these distribution functions |J. 

The other approach we have considered is an algebraic technique. In this approach, 
one derives two similar self-consistent sets of equations for the (destruction and creation) 
distribution functions by taking the advantage of the operator nature in the ground state as 
given in Eqs. (8)-(9). These self-consistent set of equations can then be tackled by various 
methods such as iteration method. As we shall see, a most simple approximation to one 
of these two self-consistent sets of equations reproduces the full CCM results, but one can 
easily go beyond that. We shall mainly discuss the algebraic approach in the followings. 

III. DISTRIBUTION FUNCTIONS AND THEIR SELF-CONSISTENT 

EQUATIONS 

We first introduce the so-called bare distribution functions as expectation value of the 
configuration operators, namely 

gi gi=(C!), (12) 

where the expectation value is defined in the usual sense as in Eq. (11). In general these bare 
distribution functions are nontrivial functions of {Fj, Fj}. Multiplying by the corresponding 
coefficients, we obtain Fjgi and Fjgi which are the the usual full distribution functions useful 
in the diagrammatic approach || . 
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Direct calculation of these functions is certainly not an easy task. Fortunately, by taking 
the advantage of the properties of the operators, one can derive self-consistent sets of equa- 
tions which can then be tackled by various methods. In particular, as C\ commutes with 
S = J2i FiCj and Cj with S = J2i FiCi, one can write 

g, = i($|e%]e 5 |$) = ±<$|eVc}|*) , (13) 
gj = I($| e %Ve s |$) = l($|C 7 eV|$> , (14) 

where A = (^\^) is the normalization constant. In order to find another expression, one 
inserts the identity e~ s e s in the expression of gi as 

gi = l($|e S Cje- 5 eV|$> = -^<$|C{eV|$) , (15) 

where the similarity-transformed operator C\ can be expanded in the nested commutator 
series as 

C\ = JC\e-~ 8 = C\ + C\\ + [5, C}] + • • • , (16) 

and this series is finite as C\ is finite and S contains only the destruction operators. By 
definition, ($|Cj = 0, hence ($|C| can be expressed in a form linear in the destruction 
operators Cj and finite order polynomial in the coefficients Fj. The expectation values of 
these finite order terms is therefore linear in gj and finite order polynomial in Fj. This 
yields a linear relation between gi and {g,j}. Hence we write 

gi = G({~ gj },{Fj}), (17) 

where G is a function linear in gj and finite order polynomial in Fj. 

In a similar fashion we write, by inserting identity e s e~ s in the expression for gi, 

gj = i($|eVCj|$> , (18) 

with the usual commutation series 

Cj = e~ s C ie s = C I + ±[C,, S] + 1[C7 7 , S], S] + ■ ■ ■ , (19) 
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and we obtain 



§i = G({gj},{Fj}), 



(20) 



where G is the same function as in Eq. (17) but now linear in {gj} and finite order polynomial 
in {Fj}. As function G is the same in Eqs. (17) and (20), only one calculation is necessary. 

Eqs. (17) and (20) provide two self-consistent sets of equations for gi and gj in terms 
of the correlation coefficients {Fj,Fj}. We note that for a particular gj its equation in 
general contains a higher-order set {gj}, and vice versa, even for a truncated coefficient 
set {Fi,Fr}. Therefore, in order to make any practical calculation, one has to make two 
approximations, a truncation on the number of coefficients {Fj,Fj}, and a truncation on 
the number of bare correlation functions {gi,gi}. After these two truncations, one should 
be able to solve the self-consistent set of equations to obtain gj and gj in terms of Fj and Fj. 
This is contrary to the CCM where one needs only one truncation (in {F I: Fj}). As we shall 
see, for a similar truncation in the coefficients {Fj, Fj}, one of our lowest order truncations 
in the bare correlation function {<?/,<?/} will reproduce the full CCM results. However, it is 
a simple step to go beyond this approximation by including some higher-order distribution 
functions gi which has proved to be essential to obtain the consistent long-range behaviors 
of the spin correlation functions (and the low-lying excitation energies) as we shall see in 
our spin model calculation. 

Since Hamiltonian usually contains terms involving both creation operators Cj and de- 
struction operators Cj, using the similar argument to the distribution discussed above, it is 
not difficult to obtain that the expectation value of a general Hamiltonian can be expressed 
as a function linear in gj and gj and a finite order polynomial in Fi or Fj, 



This expression is not unique; using Eqs. (17) or (20), one can express the Hamiltonian 
expectation as a function linear in gj and finite-order polynomial in Fj only, or as a function 
linear in gj and finite-order polynomial in Fj only, 



(H) = n({ 9l }, {~ 9I }, {Fj}) = H({~ 9I }, { 9I }, {Fj}) . 



(21) 
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(H)=H'({g I },{F I }) = H'({g I },{Fi}). (22) 

This expression of the Hamiltonian expectation is useful when we compare with the tradi- 
tional CCM. Solutions of Eqs. (17) and (20) can be substituted into these equations and we 
obtain (H) as a function of {Fi,Fi}. Variational calculation in Eq. (10) can then be carried 
out. In the following, we consider a simple application to the spin model as a demonstration. 

IV. TWO-SPIN FLIP APPROXIMATION IN A SPIN MODEL 

As a demonstration, we consider a simple truncation approximation in which the corre- 
lation operators S and S retain only the two-spin flip correlations as 

|^ 2 ) = e 52 |$), (^ 2 | = ($| e ^, (23) 

where 

fc = E/*4r' ^ = ^^4f • (24) 

Using the usual angular momentum commutations 

[sf, 4] = ±s±8 lv , [st, S; 7] = 2sf8 lv , (25) 

and the Neel state eigenequations, sf|$) = s|<3>), sj|<I>) = — s|$), it is a straightforward 
calculation to derive expectation value of various operators with respect to the states of 
Eqs. (23). In this approximation, for example, the order parameter is derived as 

M' = (8f) = 8 -'En r , (26) 

r 

where n r is the full one-body distribution function given by 

n ij = fijQij = fij ^7 ? (27) 

and we have taken the advantage of translational invariance by writing = n r with j = i+r; 
the usual two-spin correlation function is given by 
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where G^^f is the full two-body distribution function 

Gijj'ji = fijfi'j'gij,i'j' — fijfi'j' (2 s) 2 ' ' 

and finally, the expectation value of Eq. (1) is then given by 

{Hij) = -As 2 + s (ga + gij + A ]T n Vi + A ]T ™if) " A (E G H',i?i + n *i) ■ ( 30 ) 

V f i'f 

As can be seen, these physical quantities involve up to two-body distribution functions. 

The self-consistent set of equations for the bare distribution functions are derived as 
described in Sec. III. In particular, the equation for the one-body function is 

9iiji ~ fhji "I" E/ fijifhjdij 
ij 

2 

~2~fit 7i (^""1 fijidiji + E fhj9hj) 
i j 

~ >r ~2^ E fijifhjfiij'9iij,ij' + ^ E fijifi'jifhj9iji,i'j 
ijf H'j 

2 2 4 

+ /0 n2 E fijifi'jifhjfhj'9iiji,ij,i'j' ■ 

(31) 

The equation for the two-body function will contain up to twelve-body functions, 

etc. The hermitian conjugate of these equations are the self-consistent set of equations for 
9iji9ij,Vf 'j etc. Clearly, we need to make further truncations for any practical calculation. 

Consider a simple truncation in which we retain only the first two terms in Eq. (31), 
noticing that all other terms are higher-order in terms of l/2s expansion, 

9hji ~ fhji + E fijifhj9ij j (32) 

and similar equation for g^. Using the Fourier transformation technique and translational 
symmetry, it is easy to solve the two equations to obtain 

f h ~ f k /qq\ 

9k ~ i_~7~T ' 9k ~ i - f f ' ^ ' 
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where gu and fk are Fourier transformations of g^ and fy, etc. To the same order in l/2s, 
the ground-state energy, Eq. (30), is 

e = 4? E<#M+p> = - As " + + ft) + 2As ]Tn r , (34) 
where z is the number of nearest-neighbor sites and 

Elkfk ~ \ ^ 7fc/fc \ ^ \ ^ fkfk fnr\ 

- — 0i = 2^ — rr» L n ^L; — rr> ( 35 ) 

k 1 — Jkjk k 1 — Jkjk r k 1 — Jkjk 

with 

lk = -J2^' rp ■ ( 36 ) 

The variational equations, = Jj- = 0, reduce to a quadratic equation for f k and f k . The 
physical solution to these equations is 

fk = fk = ^(-1 + Jl - 7 fc 2 /A 2 ) . (37) 
7fc 

The ground-state energy and order parameter are then obtained as 



e = -As 2 + S A £(-1 + ^1 - 7fe 7 A 2 ) (38) 

and 

^ = ^~^E( / 1 2/A2 ~ 1 )- ( 39 ) 
2 k JI-tI/A 2 



It is not difficult to include the contribution of higher-order many-body distribution 
functions within our variational formalism. Consider the two-spin correlation function of 
Eq. (28), which contain the important full two-body distribution function. The bare two- 
body distribution functions gij t i>j> and gi^vy can be calculated through their self-consistent 
set of equations by keeping the same order terms in the (I /2s) expansion as we have done 
for the one-body distribution function in Eqs. (31)-(32). Without going into details of 
derivation, to the same approximation, we obtain the following results 

9ij,i'j' ~ 9ij9i'j' + 9ij'9i'j (40) 
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which, in fact, is usually referred as random-phase approximation. A simpler way to obtain 
the same results for the two-body functions is to employ the following sequential equation 
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~gj^9ij - 9ij,i'j> ~ 9ij9i'j' , (41) 

to gij and g~ij in Eq. (32) and its hermitian conjugate. 

Hence, the normalized two-spin correlation function becomes 



s„-s 



i a i+r 



)-(sl){sl +T ) = -g r g r = ~g I r . (42) 



In fact, our above results of the ground-state energy, order parameter and the correlation 
function are the same as that of the spin- wave theory ||1C|11||. In particular, the long-range 



behavior of the correlation function c r oc 1/r 2 as r — > oo for a square lattice system at 
A = 1 can not be obtained without the contribution of the two-body distribution function. 
Using iteration method for solving the equations of g^ and g^, it is straightforward to 
include higher-order contribution for other physical quantities. But we refrain ourselves 
from more detailed calculation in this article as our main purpose here is to introduce the 
new variational formalism and its comparison with the traditional CCM. 



V. COMPARISON WITH THE CCM 

In order to make a more detailed comparison with the traditional CCM, we first sum- 
marize our variational extension. We apply the HHC representation of Eqs. (8)-(9) for 
the ground state wave function of a quantum many-body system to both the ket and bra 
states with two independent, hermitian conjugate correlation coefficients {Fj, Fj} which are 
determined by the variational equations as 

d(H) _ d(H) _ 
dFj dFj 

The difficult task of expressing the Hamiltonian expectation (H) = 7i({Fj, i 7 }}) can be done 
by introducing the bare distribution functions gi = (Cj) and gi = (Ci) and solving their 
self-consistency equations 
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9I = G({~gj}, {Fj}), ~ 9l = G({ 9j }, {Fj}) . 
Using the expression of Eq. (22) for the Hamiltonian expectation 

(H)=n f ({g I },{F I }) 7 

we see that the traditional CCM is equivalent to the linear approximation in the self- 
consistensy equation for gj, namely 

9i**Fi, (43) 

for all possible values of /; the Hamiltonian expectation is then reduced to a simple form 

(^^HccudFrhiFj}) , (44) 

where function 7Yccm is linear in Fi and finite order in Fi. In our spin model calculation of 
Sec. IV, within the similar truncation involving only up to two-spin flip correlations (the so- 
called SUB2 approximation), the corresponding CCM calculation is to ignore all two-body 
and higher-order many-body distribution functions in Eq. (31). The two-spin correlation 
function thus calculated has unphysical behaviors as discussed in Ref. 8. 

In Arponen's extension of the CCM ||, while keeping the ket state in the traditional 
CCM form, the bra state, also not manifestly hermitian to ket state, is parametrized by a 
nonlinear factor as 

(|r| = (ty\e § e~ s , (45) 

where correlation operators S and S are give as in Eqs. (8)- (9). Within our variational 
formalism, as can be demonstrated, this representation of the bra state is equivalent to 
applying the random-phase approximation to the gj equation to obtain gj as finite-order 
polynormial of {fj}, 

9i « 9i({fj}) ■ (46) 

A detailed calculation in our spin model revealed that in the similar SUB2 truncation in 
Arponen's approach, the random-phase results of Eq. (40) for the two-body distribution 
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function is reproduced. While it is clearly an improvement over the traditional CCM, Arpo- 
nen's extension is also known to be poor for the strongly correlated systems (e.g., quantum 
Helium-4 fluid)even when high-order many-body contributions are considered. This may be 
related to the random-phase approximation of Eq. (46). 

In conclusion, our variational extension of the CCM provides a general many-body the- 
ory in which the traditional CCM represents a simple linear approximation and Arponen's 
extension represents a random-phase approximation. The traditional CCM is well known 
to be efficient in obtaining accurate ground state energy for a finite system with a large en- 



ergy gap separating the ground and excited states [jl^l ; but, as pointed out earlier, it is poor 
when long-range correlations in the system are important; and Arponen's extension provides 
a remedy in producing the physical long range behaviors. But both these two methods are 
known to be poor in dealing with strongly correlated systems which demand correct descrip- 
tion at short-range, in most cases, strongly repulsive. We believe the variational formalism 
represented here may provide an effective approach to the strongly correlated systems. Fur- 
thermore, its strong overlap with other well known many-body theories such as the method 
of correlated basis functionals can provide useful clues in making suitable approximations 
in practical calculations. In this regard, we are encouraged to read a recent preprint [D^] in 



which a fully variational approach has been employed to study several weakly interacting 
boson systems using the so-called independent pair correlation functions, which, in fact, 
is the corresponding SUB2 truncation of the Eq. (23) but writing in the real space. We 
like to point out that a clear advantage of the algebraic approach presented here is that it 
is straightfowardly extendible to include higher-order many-body correlations beyond the 
SUB2 level and as well as to other systems as electrons. 
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